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ABSTRACT 

A finite element code for analysis of inviscld subsonic flows over 
arbitrary non-lifting planar or axisymmetrlc bodies is described. The code 
solves a novel primitive variable formulation of the coupled irrotationallty 
and compressible continuity equations. Results for flows over a cylinder, a 
sphere, and a NACA 0012 airfoil verify the code. Computed subcrltical flows 
over an axisymmetrlc boattailed afterbody compare well with finite difference 
results and experimental data. Iterative coupling with an Integral turbulent 
boundary layer code shows strong viscous effects on the inviscld flow. Im- 
provements in code efficiency and extensions to transonic flows are discussed. 


INTRODUCTION 

Modern tactical aircraft usually have one or two jet engines within the 
fuselage, with the jet exhausts exiting through axisymmetrlc nozzles at the 
rear. Integration of the afterbody and nozzle, or boattail region (fig. 1), 
with the fuselage can strongly effect aircraft performance W . Indeed, the 
boattail drag can be as much as 20% of the total aircraft drag at transonic 
speeds [aJ . 

Several flow phenomena contribute to the boattail drag. The boattail is 
immersed in the thick turbulent boundary layer that develops over the fuse- 
lage. The flow undergoes a drag producing expansion onto the boattail followed 
by a recompression onto the jet wake, with possible shocks and separation add- 
ing to the drag. Finally, the high velocity jet displaces and entrains the 
external flow, thereby modifying the boattail pressure distribution. If the 
jet entrainment is neglected, the jet may be modeled as a solid body or sting. 

This paper describes a new finite element (FE) method for solving the 
primitive variable equations of inviscld subsonic flow about arbitrary 
axisymmetric or 2-D symmetric bodies. Subcrltical flow over a boattail 
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model, including viscous interaction using an integral boundary layer method, 
served as the primary test of the method. 


Many computational methods for the boattail problem have been published. 
Here they are grouped by the method used for the inviscid flow. An integral 
equation method for Incompressible inviscid flows was coupled with an integral 
Ir.'t. idary layer method in [2J. The classical axlsymmetric transonic full poten- 
• 1 flnlt’j difference (FD) boattail solution is given in References 

[4 and 5j built on by adding boundary layer and exhaust plume effects. 
Similar but Independent methods were used in ^6 and 7J. Reference M is an 
extension of H to three dimensions and angle of attack. The axlsymmetric 
Navier-Stokes equations were solved by FD methods in ^9 and 10^ and by a FE 
method in Ql]. A 3-D Navier-Stokes FD solution for supersonic flows over 
ogive cylinder-boattails was presen^d in ^12J. Additionally, an experimental 
study of boattail drag is given in [l3J. To our knowledge, reference [llj is 
the only other FE solution of the boattail problem. 


Many applications of FE methods to other inviscid flow problems have been 
published recently. Incompressible potential flows were considered in ,[i«] 
and (15]. Subcrltical potential flows were solved in 2-D in QI6 and 1?3 and 
in 3-D in [l8^ . References Q9 through 22 \ solved 2-D transonic potential 
flows and [53^ included 3-D transonic flows. Primitive variable FE methods 
are described in [24^ for subcrltical flows and in £25j| for transonic flows. 
The method described in [^25j is somewhat similar to the one described here. 


In the present work a FE method was used to solve a novel formulation of 
inviscid flow. The equations solved here, referred to as the nonhomogeneous 
Cauchy-Riemann equations in [26] , consist of a compressible continuity equa- 
tion and the irrotationality condition. Both are first order equations for 
the velocity components. Steady first order equations are difficult to solve 
by FD techniques but simplify the FE formulation by allowing first order in- 
terpolation functions. The form of the compressible continuity equation used 
reduces to the incompressible equation as M^d -*■ 0. Compressible flows are 
solved iteratively, starting with the Incompressible solution. The iterative 
scheme converges rapidly for subsonic flows. 


Flows over simple geometries have been computed and are presented to 
verify the method. Boattail solutions are compared to a surface source solu- 
tion for incompressible flow and to a FD solution [6] for compressible flow. 
An Integral boundary layer method [[27] was coupled iteratively with the 
inviscid calculation. Viscous boattail results are compared to experimental 
data from [l3]. Finally, methods for improving the coraputat 5 onal efficiency 
of the scheme and for extending it to tran.sonic flows are disc ussed. 


EQUATIONS AND BOUNDARY CONDITIONS 

As derived In [28], the axlsymmetric inviscid continuity equation for 
isentropic flow is: 
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v^ere all velocities are normalized by the free-stream velocity and all 
lengths by an arbitrary length scale. The axisynmetric terras are multiplied 
by a switching integer j, such that 1*1 gives the axlsymmetric equation 
and J ■■ 0 gives the two-dimensional equ^^.tion. 

This equation is solved simultaneously with the normalized irrotation- 
ality condition: 


- iX . Q (2) 

3y 3x 

The normalized boundary conditions are the free-stream conditions: 

As x,y -+■ u -► 1, V -► 0 (3) 

and the tangency condition at a solid surface: 

^nb “ vjjcos 0 J, - u^sin 65 “ 0 (4) 

where uj) and V], are the x and y velocity components on a surface at 
angle 6 ^^ to the x-axis. Equations (1) and (2) are equivalent to the non- 
conservative full potential equation. Note that (2) is linear and (1) is of 
the form: 

linear continuity equation - non-linear terms ■■ 0 

so that as 0 the equation set becomes linear. 


FINITE ELEMENT FORMULATION 

The dependent variables u and v, and for axlsymmetric flow the inde- 
pendent variable y, are approximated within arbitrary quadrilateral Fiji's by 
interpolations of the form: 

u = %Ufj V = y = Oj,yj, (5) 

where are bilinear interpolation functions dependent on element geometry, 

Uni v^i and yjq are values of u, v, and y at node N, and the summation 
convention is implied. 

Substitution of the approximations (5) into the equations of motion (1) 
and (2) results in residual errors which are minimized by the Galerkln method. 
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a subclafis of the method of weighted residuals. The Galerkln formulation of 
(2) results la: 
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which may be written as: 
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Similar application of the Galerkln method to (1) gives the system of equa- 
tions: 


C - ^^2 [F(u"-1,v'»- 1), (D + jE) - M£[H(un-l,vn-l)' 


+ G(un-l,vn-l)] 


+ jl ( u ""^, v **~^)3 


A B 

whare subscripts N and M have been omitted for clarity. 


(7) 


Equation (7) represents eight algebraic equations In eight unknowns - 
the eight velocity components Ju^.v^T at the four nodes M of an element. 

Terms A-1 are (4x4) coefficient matrlcles given In the appendix and are 
evaluated numerically by two point Gaussian quadrature. Superscripts ( )” 
refer to Iteration number. Equations (7) may be written for each element of 
a flow field, then assembled Into a global matrix equation. It can be shown 
that the assembled equations are second order accurate. 


NUMBERICAL TREATMENT OF BOUNDARY CONDITIONS 

Inflow and free stream conditions on u and v (3) are applied at a 
finite distance from the body by direct substitution Into the global matrix 
equation. Elimination of these terms would have been computationally more 
efficient but was not done. Along symmetry axes and the outflow boundary the 
V velocity component was set to zero and the u component was solved far 
from the flow equations. Resulting values of u along the outlet were within 
a fraction of a percent of the free-stream values. However, when the u com- 
ponents at the outlet were fixed at exactly their free stream values, the 
global matrix equation became singular. 

Flow tangency conditions (4) were applied using a LaGrange multiplier 
technique discussed In This technique introduces an auxiliary con- 

straint equation with a non-physical unknown at each body node. Because of 
the extra unknowns, the technique is Inefficient and would not be recommended 
for future work. 
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ITERATIVE SOLUTION METHOD 

Terns F(u,v) - I(u,v) in (7) are non-linear functions of the nodal 
unknowns and are evaluated interatively. For the first iteration (n * 1) 
the terms F - I arc set to 7.ern and the resuitinK set is solved using a 
Guass elimination scheme for banded matrices, giving the incompressible 
solution. In subsequent Iterations the non-linear terms are evaluated from 
the previous values of u and v. This scheme Is equivalent to lagging the 
density calculation, and is unconditionally stable for subsonic flows* Con- 
vergence is extremly rapid, with all - u® and v®*^ - v® less than 

10"^ in five to seven iterations. 


VISCOUS-INVISCID INTERACTION 

The FE Inviscid flow code was coupled Iteratively with a SaeeanrCresc 
[27] type integral boundary layer code by the classical method of augmentation 
of the body by the displacement thickness. Compressible inviscid calculations 
were allowed to fully converge before the viscous corrections were made. It 
was shown in that considerable computational time can be saved by using 
partially converged inviscid solutions to update the boundary layer, but this 
was not attempted. 


COMPUTATIONAL GRID 

Sheared grids, as illustrated at the tops of figures 4 and 6, were used 
for each FE solution mentioned below. The bandwidth of the global matrix that 
must be stored and solved is proportional to the number of nodes along a vert- 
ical grid line; so the vertical extent of the grid was kept small. Solution 
"wiggles" were often noted on highly elongated elements. These could presum- 
ably be avoided by using a less distorted grid such as the wrap around grids 
described in [303* In this paper, solution wiggles have been minimized by 
plotting flow quantities evaluated at element midpoints using interpolation 
equations (5). 


VERIFICATION RESULTS 

Figures 2 and 3 show calculated pressure coefficients on a sphere and a 
cylinder, respectively. The symmetric solutions are shown only over the 
second quadrant, although both solutions were computed over the upper half 
plane. The same coarse 9x53 node grid (17 points on the body) was used for 
both calculations. Incompressible FE solutions agree well with the analytic 
solutions. Compressible results at ■ 0.3 for the sphere and Mod ” 0.38 
for the cylinder, both near their critical Mach numbers, agree reasonably well 
with those predicted by a Gothert rule compressibility transformation of the 
analytic incompressible solution. 

Flow over a non-lifting NACA 0012 airfoil at Mo> * 0.7 was computed on 
the 18x47 node grid shown at the top of figure 4. Nineteen points were dis- 
tributed over the body at the locations tabulated in Abbot and VonDoenhoff 
^31j. The lower half of figure 4 shows the computed Mach number contours at 
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VL ■ 0.7. Computed incompressible and compressible pressure coefficients are 
compared with experimental data from [31 and 323 respectively in figure S» 
with good agreement in both cases. 


BOATiAIL RESULTS 

Computed pressure distributions for a boattail model are compared with 
experimental data from [l33» which also gives details uf the geometry. 

Briefly, the boattail model, which is sketched in figure 7, starts with a 
10° half angle conical nose followed by a long cylindrical centerbody. The 
boattail region, which makes up about nine percent of the model (excluding 
the sting) consists of a spherical section that smoothly joins the centerbody 
to a 15° conically necked-down se^'. £ion. The model is mounted on a long cylin- 
drical sting representing a propulsive jet. 

Figure 6 shows a 10x102 node grid (85 points on the body, nine over the 
boattail region) used for the boattail computations. The lower half of the 
figure shows computed Mach number contours at Mn ** 0.8. 

Incompressible results for the boattail are compared with a surface 
source solution in figure 7. Agreement is good over the entire body except 
near the nose. A good incompressible FE solution is important since it acts 
as a starting point of the compressible Iteration. 

A subcritical solution over the boattail region, both with and without 
viscous interaction, is compared with experimental data in figure 8. The 
tree stream Mach number is Moo “ 0.7 and the Reynolds number based on body 
length excluding sting is Rcl * l,3A13xl07, The viscous solution underpre- 
dicts the initial expansion but agrees well with the data in the recompres- 
sion over the boattail and onto the sting. The figure illustrates the strong 
effect that a thick boundary layer can have on an invlscid pressure distribu- 
tion. At the bottom of the figure the actual body is shown with the equivalent 
displacement body, illustrating the smoothing effect the boundary layer has on 
the boattail-sting juncture. 

A small trick was needed to compute the viscous flows shown in figures 8 
and 9. In the first iteration of the invlscid-viscous interaction, the bound- 
ary layer code invarably predicted separation in response to the strong 
invlscid pressure gradient near the end of the boattail. However, the con- 
verged interaction shows no separation. To avoid separation at the first 
Iteration, the invlscid solution was first converged at M» “ 0.2. The bound- 
ary layer did not separate in response to this pressure field, so an equivalent 
displacement body could be determined. Subsequent invlscid solutions at the 
desired Mach number but over the current displacement body were enough smoothed 
that the boundary layers did not separate, and the invlscid-viscous interaction 
could be converged. 

Another boattail solution at ^ 0.8 and Re^ “ 1.4181x10^ is shown in 
figure 9. (Mach number contours ft : the invlscid solution are shown in fig. 6.) 
The invlscid solution is Interesting in tiat it shows supersonic velocities at 
four nodes. The FE formulation is purely elliptic and contains no upwind bias 
in supercritical regions. Without special treatment FD relaxation schemes 
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typically diverge for supercritical flows. In fact* the FE solution also 
diverges at slightly higher It is thought that the direct solution of 

the FE equations over the «itire flowfield is responsible for the method's 
ability to converge with slightly supercritical flow. 

Also shown in figure 9 is Chow, Bober, and Anderson's FD non-conserva- 
tive transonic potential solution for the same problem. The FD solution, 
which includes upwind differencing in supersonic regions, predicts a strong 
shock recompression on the boattail. This shock is completely eliminated 
when the viscous interaction is considered. The FE and FD codes use the 
same boundary layer code, but the FD solution shows better agreement with 
the viscous data. This may be due to differences in treatment of the free- 
stream boundary, which is mapped to Infinity in the FD code but Is placed 
about 10 body radii away in the FE code. 


COMPUTER REQUIREMENTS 

Both the FE and FD codes were run on an IBM 370/158 computer. CPU times 
for invlscld solutions like those shown in figure 9 were 760 seconds for the 
FE solution and 260 seconds for the FD solution on similar grids. The FE 
code has since been cleaned up and run on an IBM 3033. The invlscld solution 
in figure 9 took about 34 seconds and the interaction solution took about 
170 seconds. The FE solution required over one million words of storage, 
about 13 times that required by the FD code, principally for storage of the 
global coefficient matrices. Storage requirements are not much of a problem 
on a virtual storage computer, as page faults can be minimized by accessing 
the global matrix only columnwise. 


DISCUSSION AND CONCLUSIONS 

A FE solution of the first order equations of invlscld, Irrotational 
compressible flow has been developed, verified, and applied to the analysis 
of subsonic flow over a boattail. Reasonable agreement with experimental 
boattail data was obtained when the Invlscld code was coupled with an Integral 
boundary layer code. Two-dimensional Inviscld-vlscous interaction calcula- 
tions were made quickly, and thus appear to be useful tools for screening 
boattailed afterbody designs prior to more extensive calculations or wind 
tunnel testing. 

Certain aspects of this FE algorithm are appealing. First, FE methods 
are able to fit arbitrary boundaries, with the mapping to computational space 
contained In the element formulation rather than In transformation of the flow 
equations. Incidentally, the Interpolation functions for the flow variables 
are also useful when flow variables are needed at off-node points, for example 
when calculating a detailed pressure distribution for boundary layer coupling. 
Secondly, the primitive variable formulation allows use of bilinear Interpola- 
tion functions, guarantees a continuous velocity and density solution, and 
eliminates the need to differentiate a potential function. However, solving 
two equations per node lengthens solution time. Finally, the ability of the 
method to model the first order primitive variable equations even at sonic 
velocities and without convergence problems suggests that direct solvers may 
be useful for other two-dimensional problems. 
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The FE code was slo%^r than an existing FD coda» hut CPU times for the 
FE code could be reduced by code refinement. Time spent in Gaussian solution 
of the global matrix equations could be reduced by vectorization for a vector 
computer. 

The artificial compressibility technique developed by Hafez [l9] could 
be used to extend the present method to transonic flows. Here the continuity 
equation would be written in conservative form. The density calculation is 
lagged, and upwinded at supersonic nodes. Since the conservative form of the 
continuity equation has many fewer terms than (1). the artificial compressi* 
bllity technique should also improve the computational efficiency of the 
method. 
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APPENDIX 

Terms AnM INM in equation (7) are (4x4) coefficient matrices result- 
ing from the Galerkln FE formulation of the irrotationality and continuity 
equations, and are given by: 

anm • y % ^ 

A 

% y^dydx 
A 

CnM • “ ®NM 
®NM * 

®NM ° / %%dydx 

A 

A 

“ 2 J' «Ny-^dydx 

HnM • y + (.QlVl){|% y-’dydx 

A 
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Inm ■ f “ ^%tlMdydx 

A 

K,L«M,N - 1,2,3. 4 


Terms A)^ - are functions of known Interpolation functions Q}]. 

Terms also Include the nodal unknowns ui^ and vj^ and must be 

evaluated Iteratively. All Integrals were evaluated numerically In the 
present work; however, computer time could probably be saved by evaluating 
some of the Integrals analytically. 


REFERENCES 

1. Hunt, B. L., Surber, L. E. . and Grant, G. K. , "Propulsion-System Inte- 

gration for Tactical Aircraft," Astronautics and Aeronautics , Vol. 18, 
June 1980, pp. 62-68. 

2. Nakayama, A., Patel, V. C., and Landweber, L., "Flow Interaction Near 

the Tail of a Body of Revolution, I - Flow Exterior to Boundary Layer 
and Wake," Journal of Fluids EnRineering . Vol. 98, Sept. 1976, 
pp. 531-537; "II - Iterative Solution for Flow Within and Exterior to 
Boundary Layer and Wake," Journal of Fluids Engineering . Vol, 98, 

Sept. 1976, pp. 531-546. 

3. South, J. C. , Jr. and Jameson, A., "Relaxation Solutions for Invlscld 

Axlsymmetrlc Transonic Flow Over Blunt or Pointed Bodies," AIAA Compu- 
tational Fluid Dynamics Proceedings, Palm Springs, Cal., July 9, 1973. 

4. Cosner, R. R. and Bower, W. W. , "A Patched Solution of the Transonic 

Flowflelds About an Axlsymmetrlc Boattall," AIAA Paper 77-227, Jan. 1977. 

5. Kuhn, G. D. , "Computer Program for Calculation of Separated Turbulent 

Flows on Axlsymmetrlc Afterbodies," Arnold Engineering Development Center 
Report AEDC-TR-77-72, June 1977. 

6. Chow, W. L., Bober, L. J., and Anderson, B. H. , "Numerical Calculation of 

Transonic Boattall Flow," NASA TN D-7964, June 1975. 

7. Wilmoth, R. G., "Computation of Transonic Boattall Flow with Separation," 

NASA TP-1070, 1977. 

8. Nakayama, A. and Chow, W. L. , "Calculations of Transonic Boattall Flow at 

Small Angle of Attack," University of Illinois, Champa Ign-Urbana, 111. 
Report ME-TN-395-6, April 1979 (NASA CR-158471). 

9. Holst, T. L. , "Numerical Solution of Axisymmerrlc Boattall Flow Fields 

with Plume Simulators," AIAA Paper 77-224, Jan. 1977. 

10. Cosner, R. R. , "Fast Navier-Stokes Solution of Transonic Flowflelds About 
Axlsymmetrlc Afterbodies," AIAA Paper 80-0193, Jan, 1980. 


10 


11. Cooke. C. H. and Blanchard. D. K.. "A Shock Capturing implication of the 

Finite Element Method*" International Journal for Numerical Methods in 
Engineering . Vol. 14. 1979. pp. 271-286. 

12. Schiff, L. B. and Sturek. W. B.. "Numerical Simulation of Steady Super- 

sonic Flow Over an Oglve-Cyllnder-Boattall Body," AIAA Paper 80-0066. 
Jan* 1980* 

13. Shrewsbury, C. D*. "Effect of Boattall Juncture Shape on Pressure Drag 

Coefficient of Isolated Afterbodies," NASA TM X-1S17, 1968. 

14. Baskharone, E. and Hamed, A., "A New Approach to Cascade Flow Analysis 

Using the Finite Element Method," AIAA Paper 80-0389, Jan. 1980. 

15. Keck, H. and Haas, W. . "Vorticlty and Blade Circulation Modelling In the 

Calculation of Quasi-Three-Dimensional Cascade Flows with Finite Ele- 
ments," Third International Conference on Finite Elements in Flow 
Problems, Banff, Alberta, Canada, June 10-13, 1980. 

16. Chima, R. V., "Finite Element Analysis of Inviscid Subsonic Flow with 

Viscous Interaction and Inverse Applications," Ph.D. Dissertation, 
University of Akron, 1979. (Available from University Microfilms.) 

17. Habashi, W. G., "Finite Element Approach to Compressor Blade-to-Blade 

Cascade Analysis, " AIAA Journal . Vol* 27, July 1979, pp. 693-698. 

18. Laskaris, T. E., "Finite Element Analysis of Three-Dimensional Potential 

Flow in Turboraachines," AIAA Journal . Vol. 16, July 1978, pp. 717-722. 

19. Hafez, M. S. and Murman, E. M. , "Artificial Compressibility Methods for 

Numerical Solution of Transonic Full Potential Equation," AIAA Paper 
78-1148, July 1978. 

20. Ecer, A. and Akay, H. U., "Investigation of Transonic Flow in a Cascade 

Using an Adaptive Mesh," AIAA Paper 80-1430, July 1980. 

21. Deconlnck, H. and Hlrsch, C. H., "A Finite Element Method Solving the 

Full Potential Equation with Boundary Layer Interaction in Transonic 
Cascade Flow," AIAA Paper 79-0132, Jan, 1979. 

22. Brlsteau, M. 0. et al., "Transonic Flow Simulations by Finite Elements 

and Least Square Methods," Third International Conference on Finite 
Elements In Flow Problems, Banff, Alberta, Canada, June 10-13, 1980. 

23. Eberle, A., "Finite Element Methods for the Solution of Full Potential 

Equation In Transonic Steady and Unsteady Flow," Third International 
Conference on Finite Elements In Flow Problems, Banff, Alberta, Canada, 
June 10-13, 1980. 

24. Fletcher, C. A. "A Primitive Variable Finite Element Formulation for 

Inviscid, Compressible Flow," Journal of Computational Physics, Vol. 33, 
1979, pp. 301-312. 


11 


25* QuCCoe» J. J.» Gulu-Roux> J.» and Laminlc, J., "Finlt* ElananC Calcula- 
Cion of Steady Transonic Flow in Nozales Using Prinary Varlablss/' 

7th Intsmatlonal Confsrsncs on Numarlcal Hsthods In Fluid Dynsnics* 
University of Stanford, Stanford, Calif., June 23-27, 1960. 

26. Lomax, H. and Martin, E. D., **Fast Direct Numerical Solution of the Non- 

homogeneous Cauchy-Relmann Equations," Journal of Computatlorsal Physics, 
Vol. 15, 1974, pp. 55-80. 

27. Sasman, P. K. and Crescl, R. J., "Conq>resslble Turbulent Boundary Layer 

with Pressure Gradient and Heat Transfer," AIAA Journal . Vol. 4, Jan. 
1966, pp. 19-25. 

28. Liepmann, H. W. and Roshko, A., Elements o f Gas Dynamics . Wiley, New York, 

1957. 


29. Chung, T. J., Finite Element Analysis In Fluid Dynamics . McGraw-Hill, 

New York, 1978. 

30. Sorenson, R. L. , "A Computer Program to Generate Two-Dimensional Grids 

About Airfoils and Other Shapes by the Use of Poisson's Equation," 
NASA TM 81198, May 1980. 

31. Abbott, I. H. and von Doenhoff, A. E., Theory of Wing Secti ons , Dover, 

New York, 1959. 


32. Experimental Data Base for Computer Program Assessment," AG4RD AR-138, 
May 1979. 







Figured. - Top-. NACA 0012 airfoil grid. Bottom; Mach number 
contours. Mi,o-a7. 



Figure 6. -To(k Boattailgrid. Bottom: Mach number contours, ^*0.8. 







